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Abstract In the physically non viscous fluid dynam- 
ics, "shock capturing" methods adopt either an artificial 
viscosity contribution or an appropriate Riemann solver 
algorithm. These techniques are necessary to solve the 
strictly hyperbolic Euler equations if flow discontinuities 
(the Riemann problem) must be solved. A necessary dis- 
sipation is normally used in such cases. An explicit artifi- 
cial viscosity contribution is normally adopted to smooth 
out spurious heating and to treat transport phenomena. 
Such a treatment of inviscid flows is also widely adopted 
in the Smooth Particle Hydrodynamics (SPH) finite vol- 
ume free Lagrangian scheme. In other cases, the intrinsic 
dissipation of Godunov - type methods is implicitly use- 
ful. Instead "shock tracking" methods normally use the 
Rankine - Hugoniot jump conditions to solve such prob- 
lem. A simple, effective solution of the Riemann prob- 
lem in inviscid ideal gases is here proposed, based on an 
empirical reformulation of the equation of state (EoS) 
in the Euler equations in fluid dynamics, whose limit 
for a motionless gas coincides with the classical EoS of 
ideal gases. The application of such effective solution of 
the Riemann problem excludes any dependence, in the 
transport phenomena, on particle resolution length h in 
non viscous SPH flows. Results on ID shock tube tests 
are here shown. 
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1 Introduction 



In both Lagrangian and Eulerian inviscid fluid dynam- 
ics, a small dissipation is normally useful to handle dis- 
continuities in the flow (the Riemann problem). An ar- 
tificial viscosity is introduced in SPH, as a shock cap- 
turing method, to prevent particle interpenetration and 
to smooth out spurious heating in the flow to solve the 
strictly hyperbolic system of Euler equations. The in- 
troduction of such a small dissipation, to solve the Eu- 
ler equations, is also currently adopted to produce both 
mass and angular momentum transport in SPH physi- 
cally inviscid modelling of accretion discs in astrophysics 
(Molteni et al. 1991, 1994; Lanzafame et al. 1992, 1993, 
1994, 2000, 2001; Belvedere et al. 1993; Chakrabarti & 
Molteni 1993; Meglicki et al. 1993; Murray 1996; Lan- 
zafame & Belvedere 1997, 1998; Okazaki et al. 2002). 
Efforts were accomplished in SPH to solve both the " ap- 
proximate" and the "exact" Riemann problem, either via 
a reformulation of the artificial viscosity term (Balsara 
1995; Monaghan 1997; Morris & Monaghan 1997) or via 
Godunov techniques (Yukawa et al. 1997; Inutsuka 2002; 
Parshikov & Medin 2002 Molteni & Bilello 2003) which 
also include an " intrinsic" dissipation. In the first case, a 
reformulation of the artificial viscosity could be necessary 
because, for "weak shocks" or low Mach numbers, the 
fluid becomes "too viscous" and angular momentum and 
vorticity could be non - physically transferred. A tech- 
nique to solve the Riemann problem, based on a refor- 
mulation of the EoS in the Euler equations, is here pre- 
sented, where particle SPH pressure terms are recalcu- 
lated without any artificial viscosity contribution. Since 
shock flows are non-equilibrium events of fluid dynamics, 
we pay attention to the fact that the EoS: p — (7 — l)pe 
for ideal flows cannot exactly be applied to solve the Rie- 
mann problem. In fact such equation is strictly valid only 
for equilibrium or for quasi-equilibrium thermodynamic 
state, when the velocity of propagation of perturbations 
equals the sound velocity. Successful results, based on 
some form of mathematical dissipation introduced within 
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the strictly hyperbolic system of Euler equations for ideal 
flows, are obtained due to the necessity to correct such 
congenital defect. In §2 of this paper, we briefly remind 
how does SPH method work for ideal non viscous flows. 
In the same §2 we also outline which evolution has been 
accomplished in the explicit artificial viscosity dissipa- 
tion description. Instead, in §3, we show how to refor- 
mulate the EoS, according to the Riemann problem for 
inviscid, ideal gases. Applications to ID shock tubes (Sod 
1978), to solve shocks are also shown in §4. 



where the sum is extended to all particles included 



within the interpolation domain D, 



pj/rrij is the 



number density relative to the jth particle. W(r i5 r ■, h) < 
1 is the adopted interpolation Kernel whose value is de- 
termined by the relative distance between particles i and 
j. Typically, such cubic spline Kernels W(rij,h) are in 
the form: 



1 - 



W ij = ^\ 1 I (2-r ij ) 3 ifl<r y 
otherwise 



< 1 

< 2 



(8) 



2 The Euler equations and their SPH 
formulation 

2.1 SPH and artificial viscosity for non - viscous flows 

As for Lagrangian ideal non - viscous gas hydrodynamics, 
the relevant equations (Euler equations) are: 



|+PV., = 



dv 



Vp 

p 



de p 

—— = V • v 

at p 

P=(l- l)pe 
dr 

—r= v 
dt ~ 



continuity equation(l) 
momentum equation(2) 

energy equation(3) 
perfect gas equation(4) 
kinematic equation(5) 



The most of the adopted symbols have the usual 
meaning: d/dt stands for the Lagrangian derivative, p 
is the gas density, e is the thermal energy per unit mass. 
The adiabatic index 7 has the meaning of a numerical 
parameter whose value lies in the range between 1 and 
5/3, in principle. 

The SPH method is a Lagrangian scheme that dis- 
cretizes the fluid into moving interacting and interpo- 
lating domains called " particles" . All particles move ac- 
cording to pressure and body forces. The method makes 
use of a Kernel W useful to interpolate a physical quan- 
tity A(r) related to a gas particle at position r according 
to (Monaghan 1986, 1992): 



MD = / A(r')W(r,r',h)dr' 



(6) 



D 



W(r,r',h), the interpolation Kernel, is a continuous 
function - or two connecting continuous functions whose 
derivatives are continuous even at the connecting point 
- defined in the spatial range 2h, whose limit for h — > 
is the Dirac delta distribution function. All physical 
quantities are described as extensive properties smoothly 
distributed in space and computed by interpolation at r. 
In SPH terms we write: 

N . N . 



even though also Gaussian form are adopted. In ex- 
pression (8), Vij = — Lj\/h represents the module of 
the radial distance between particles i and j in units of 
h. 

In SPH formalism, equations (2) and (3) take the 
form: 



dv.j 
dt 




dt -2^\pl + p - 



(9) 
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where = v t —Vj and mj is the mass of jth particle. 

For a better energy conservation, the total energy 
E = (e + \ v 2 ) can also be introduced in the SPH formu- 
lation: 

d 



dt 



-Ei 




+ 



PjVj 



of the energy equation: 
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= --V • (pv) ■ 
P 



(11) 



(12) 



^j- = J2m j v ij -V i W ij 
3=1 



In this scheme the continuity equation takes the form: 

N 



or, as we adopt, it can be written as: 

N 

™ . T/f/ 



Pi = £ m 3 W i 
3 = 1 



(13) 



(14) 



which identifies the natural space interpolation of 
particle densities according to equation (7). 

The pressure term includes the artificial viscosity con- 
tribution given by Monaghan (1985, 1992) and Mon- 
aghan & Lattanzio (1985), with an appropriate ther- 
mal diffusion term which reduces shock fluctuations. It 
is given by: 



3=1 



3=1 



Vij = a ^3 



(15) 
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where 



P-ij = < 



if Vu ■ r,,- < 



otherwise 



(16) 



with c s i being the sound speed of the ith particle, 
r t j = r f — r ■, £ 2 <C h 2 , a w 1 and /J « 2. These a and 
/3 parameters of the order of the unity (Lattanzio et al. 
1985) arc usually adopted to damp oscillations past high 
Mach number shock fronts developed by non-linear in- 
stabilities (Boris & Book 1973). SPH method, like other 
finite difference schemes, is far from the continuum limit. 
The linear a term is based on the viscosity of a gas. 
The quadratic (j3, Von Neumann-Richtmyer-like) artifi- 
cial viscosity term is necessary to handle strong shocks. 
Linear a and quadratic f3 artificial viscosity terms (usu- 
ally ~ 1 and sometimes, in some specific cases, < 1) are 
chosen = 1 and = 2, respectively. In the physically in- 
viscid SPH gas dynamics, angular momentum transport 
is mainly due to the artificial viscosity included in the 
pressure terms as: 
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(17) 



where terms into parentheses refer to intrinsic gas 
properties. 

In SPH conversion (eqs. 9, 10, 11, 13, 14) of mathe- 
matical equations (eqs. 1 to eq. 4) there are two princi- 
ples embedded. Each SPH particle is an extended, spher- 
ically symmetric domain where any physical quantity / 
has a density profile fW(r i ,r_j,h) = fW{\r_i — rj\,h). 
Besides, the fluid quantity / at the position of each SPH 
particle could be interpreted by filtering the particle data 
for f(r) with a single windowing function whose width 
is h. So doing, fluid data are considered isotropically 
smoothed all around each particle along a length scale 
h. Therefore, according to such two concepts, the SPH 
value of the physical quantity / is both the superposition 
of the overlapping extended profiles of all particles and 
the overlapping of the closest smooth density profiles of 



2.2 Progresses on artificial viscosity formulation 

Strong shock require a = 1. However for weak shocks 
and for small Mach number flows, the fluid becomes "too 
viscous" and both angular momentum and vorticity are 
transferred unphysically. To solve such problem, several 
solutions are proposed: 

a) the formulation of a "limiter" (Balsara 1995), mul- 
tiplying the artificial viscosity terms r\i for = 0.5(/, + 



fj), where 



IV -Vi 



IV-^I + IVxi^I + IO-^/ZT 



(18) 



It reduces the unhysical spread of angular momentum 
in whirling flows up to 20 times. It is w 1 for planar 
shocks, while it increases if the tangential kinematics is 
relevant; 

b) a switch to reduce artificial viscosity (Morris & 
Monaghan 1997). In this hypothesis, for each ith particle 
the a evolves according to a decay equation: 



d(Xi 



o.i — a 



(19) 



0.1, Tj = h/c s and the source term Si = 
i,0). 



where a 
fi max(— V • v 

c) a reformulation of the artificial viscosity according 
to the Riemann problem (Monaghan 1997). Particles i 
and j are considered as the equivalent left " Z" and right 
"r" states of a given contact interface. The ID Riemann 
problem is taken into account along the line joining them. 
Being the Euler equations in conservative form: 



dt dx 

the simplest Euler technique of integration is: 



At 



where numerical fluxes (Marti et al. 1991) 



r(s h s r ) = 0.5 [/; + /;-£ ia;^.| 



(20) 



(21) 



(22) 



where the e_j are the eigenvectors of the Jacobian ma- 
trix A = df /ds and A* is an average of A for the " Z" and 
"r" states. Aco* are the "jumps" of s across the charac- 
teristics: 



s r - s, - J2 Ai0 *& 



(23) 



For ID ideal flows, the eigenvalues are v, v + c s and 
v — c s , where the two including the sound velocity are the 
velocities of propagation of sound referred to the frame 
whose fluid velocity is v. Assuming that the jump in the 
velocity across characteristics could physically be ■ 
Lij/\Lij\ an d that a signal velocity v sig corresponds to 
the above eigenvalues, |A*|Z\a;* corresponds to Vgig^jV^ ■ 
Lij/\Lij\- So doing, the artificial pressure contribution in 
the momentum equation is: 




Vij = - 



KVsig,ijVij-Lij/\Lij\ 

Pij 



(24) 



where K is an arbitrary constant w 1 and pij = 
0.5(pi + pj). 

As far as the energy equation is concerned, the SPH 
formulation in this case is: 
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dt 



N 

■E 1 



gas 



Pij 



(25) 



where e*,- 



e* and e* 



0.5fe-r J1 /|r J7 |) 2 + 6 J . 



The signal velocities for the ID Lagrangian Riemann 
problem on ideal flows are reported in (Whitehurst 1995), 
on the base of Gottlieb & Groth (1988) and Toro (1992) 
results. The pressure on the contact interface p* links 
the left and the right states. For systems with one shock 
and one rarefaction wave p* € [pi , p r ] and 



P 



c s i + c sr + (vi - v r )(-f - l)/2 



127/(7-1) 



Csl/P, 



(7-l)/2 7 



+ c sr /p 



( 7 -l)/2 7 



(26) 



In the case of two shocks, a more complicated rela- 
tion (Gottlieb & Groth 1988) also involves the velocity 
v* on the contact interface. However, for some practical 
purposes, p* = ma,x(pi,p r ) and v* — 0.5(u; +v r ). Notice 
that in the presence of two very strong rarafaction waves 
p* = 0. 

The Lagrangian two velocities of waves spawned by 
the interface are: 



A; 



vi 



Csl[l+ h - 1) % /P '- 1) ] if P */pi>l 



Vl - C S 1 



Hp*/ Pi < 1 



(27) 



and 



A r — 



V r + C, r [l + (l^M^ll ] if p * /pr > 1 



lip* /Pr < 1 



(28) 



In the Lagrangian description, being v S i g the speed 
of propagation of perturbation from i to j and vice- versa 
(I o r), 



J stg,i^tj 



>i Hi ' Lij/\Lij\ 

= -Csj-Vj-rij/lZijl 



J sig 



J stg,i—^j 



C s i + Csj Eij • LijI [Lij I 



(29) 
(30) 
(31) 
(32) 



having considered the versus going from i to j. By 
performing some numerical checks, Monaghan (1997) also 
suggested other similar algebraic expressions. 



2.3 Some cautions in using artificial viscosities 

In Shu (1992), some cautionary remarks are reported 
on the adoption of artificial dissipation. In particular, 
"Because the magnitude of the viscous term does not 
affect the net shock jump conditions, many numerical 



schemes implicitly or explicitly incorporate the trick of 
artificial viscosity for halting the ever-growing steepening 
tendency produced by nonlinear effects, thereby gaining 
the automatic insertion of shock waves wherever they 
are needed (in time-dependent calculations)." However, 
the numerical viscous term should be large enough to 
spread out shock transitions only over a few resolution 
length, making shock waves resolvable. In this sense, any 
mathematical dissipation should be considered a useful 
mathematical "trick" to get correct results. 



3 How EoS matches the Riemann problem 

The solution of the Riemann problem is obtained at 
the interparticle contact points among particles, where 
a pressure and a velocity, relative to the flow disconti- 
nuity, are computed. This is also clearly shown in In- 
utsuka (2002); Parshikov & Mcdin (2002); Molteni & 
Bilello (2993), where the new pressure p* and velocity 
v* are reintroduced in the Euler equations instead of p 
and v to obtain the new solutions compatible with in- 
viscid flow discontinuities. In SPH, we pay attention in 
particular (Inutsuka 2002; Parshikov & Medin 2002) to 
the particle pressures pi and pj, in the SPH formula- 
tion of the momentum and energy equations (8) and (9), 
whose substitution with pressures p* and p*, solutions 
of the Riemann problem, excludes any artificial viscosity 
adoption, although a dissipation, implicitly introduced, 
is necessary. Therefore, the solution of the Lagrangian 
Riemann problem, only in its pressure computation, is 
enough to interface SPH with the Riemann problem so- 
lution. 

The key concept is that the physical interpretation 
of the application of the artificial viscosity contribution 
in the pressure terms corresponds to a reformulation of 
the EoS for inviscid ideal gases, whose equation: 



P 



gas 



(7 - l)pe 



(33) 



is strictly applied in fluid dynamics when the gas 
components do not collide with each other. In the case 
of gas collisions, it modifies in: 



p* = (7 — l)pe + other. 



(34) 



The further term takes into account the velocity of 
perturbation propagation (Monaghan 1997). This veloc- 
ity equals the ideal gas sound velocity c s whenever we 
treat non collisional gases in equilibrium or in the case of 
rarefaction waves. Instead, it includes the "compression 
velocity": •r i j/\r i j\ in the case of shocks (eq. 29). In 
the first case, we write the EoS for inviscid ideal gases 
as: 



p\ 



gas 



?c 2 

7 " 



(35) 



where c s — (jp/p) 1 ^ 2 = [7(7 — l)e] 1 / 2 . Instead, in the 
second case, the new formulation for the EoS is obtained 
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squaring v si g^ , so that c 2 (l - v shock /c s ) 2 is an energy 
per unit mass in the case of shock. Hence: 



P \ 



if v. 



shock 



if v 



shock 



< o 



> 



(36) 



In the SPH scheme, being: 
- 2 

Pi 2 / 1 Vshock,t 
— C ai 1 



7 



Vshock,i 



(37) 
(38) 







otherwise. 



This formulation introduces the "shock pressure term" 
p{v 2 hock —2v s hockC s )M 1 whose linear and quadratic power 
dependence on v_ij ■ Lij/\Lij \ is analogue to both the lin- 
ear and the quadratic components of the artificial vis- 
cosity terms (15). The linear term oc c s v s h oc k is based on 
the viscosity of a gas. The quadratic term oc (c s v s hock) 
(Von Neumann - Richtmyer - like) is necessary to handle 
strong shocks. These contributions involve a dissipative 
power, whose effect corresponds to an increase of the gas 
pressure. Therefore, we adopt the formulation (eq. 21) 
as p* and p* in the SPH formulation of the momentum 
(eq. 9) and energy equations (eqs. 10 or 11): 



dv.j 
dt 

~dt 



N 

= -E' 

i N 
-IE' 




Pi 



dt 




(39) 



(40) 



(41) 



This simple reformulation, allows us to to keep the 
same Courant - Fricdrichs - Lewy condition as for the 
timcstcp computation, substituting only the sound ve- 
locity c s with c s - v shock . 



4 ID Sod shock tube tests 

In this section a comparison of analytical and our SPH - 
Riemann ID shock tube test results (Sod 1978), also with 
the initial particle configuration (time T = 0), is made. 
Notice that the so called analytical solution of the Rie- 
mann problem is obtained through iterative procedures 
left-right, applying to the discontinuity the Rankinc - 
Hugoniot "jump" solution. Figg. 1 and 2 shows results 
concerning the particle density, thermal energy per unit 
mass and velocity at the same final computational time 
(T = 100). Throughout our SPH- Riemann simulations, 
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Fig. 1 ID shock tube tests as far as both analytical (solid 
line) and our SPH-Riemann (dots) results are concerned. 
Density p, thermal energy e and velocity v are plotted at 
time T = 100. Density and thermal energy of particles ini- 
tially at rest at time T = are also reported. The initial 
velocity is zero throughout. 



the initial particle resolution length is h = 5 ■ 10~ 2 . The 
whole computational domain is built up with 2001 par- 
ticles from X = to X = 100, whose mass is different, 
according to the shock initial position. At time T = 
all particles are motionless. 7 = 5/3, while the ratios 
P1/P2 = 3 and ei/e 2 = 2 in Fig. 1, and P1/P2 — 3 and 
ei/e 2 = 1 in Fig. 2, between the two sides left-right. The 
first 5 and the last 5 particles of the ID computational 
domain, keep fixed positions and do not move. The choice 
of the final computational time is totally arbitrary, since 
the shock progresses in time. 

Results, where our SPH-Riemann solution is applied 
to SPH, arc in a good comparison with the analytical 
solution. Discrepancies involve only 4 particles at most. 
This means that, according to the cautionary remarks on 
§2.3, the physical dissipation introduced in the EoS (eqs. 
37, 38) is effective. Any artificial viscosity contribution 
is totally absent, as well as no thermodynamic proper- 
ties of neighbour particle are involved, as the application 
of Godunov - type solver (Yukawa 1997; Inutsuka 2002; 
Parshikov & Mcdin 2002; Molteni & Bilello 2003) does. 



5 Discussion and conclusions 

The comparison of ID Sod shock tube analytical results 
with SPH - Riemann ones, where a revision of the EoS 
within the Riemann problem is made, are quite success- 
ful. In particular, in our modelling, neither a parametrized 
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Fig. 2 The same as in Fig. 1. In this example, the initial dis- 
continuity does not affect the thermal energy per unit mass, 
being the gas initially isothermal, while the initial velocity is 
zero throughout. 



artificial viscosity, nor any dependence on spatial reso- 
lution length h, nor a Godunov solver, nor additional 
computational time are involved. However, to give a gen- 
eralization, we need only one general EoS and not a sep- 
aration of the EoS according to the kinematic of the 
flow. To this purpose, we can generalize the EoS: p* = 
pc 2 s (l - v shock /c s ) 2 /j as: 
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l-C 



(42) 



where C — > 1 for vr — y_ l} -r^^r^ | < 0, whilst C 
otherwise. A simple empirical formulation can be, 



C 



1 / v R 

— arccot R — 

7T V C. 



(43) 



where R 3> 1. R is a large number describing how 
much the flow description corresponds to that of an ideal 



gas. To this purpose, R ~ X/d, being A 



oc p 



-V3 the 



molecular mean free path, and d the mean linear di- 
mension of gas molecules. Further, improvements can be 
made next future. The comparison with analytical re- 
sults of ID Sod shock tubes in the same test in §4 gives 
us the same results showed in Figg. 1 and 2. The sim- 
ple EoS for inviscid ideal gases: p = (7 — l)pe cannot 
be strictly applied in shock problems by the fact that, 
solving the Euler equations, not only instabilities and 
spurious heating come out, but also that this empirical 
EoS derived by Charles, Volta Gay - Lussac and Boyle 
- Mariotte laws, should be strictly applied only cither 



in equilibrium configurations or to "quasi-static" evo- 
lutions of thermodynamic systems without any dissipa- 
tion. This is a restriction that does not match with shock 
flow dynamics, when dissipation on the shock front can- 
not be neglected. Those techniques involving Godunov 
- type schemes also introduce some dissipation mecha- 
nisms (Park & Kwon 2003). Therefore, according to these 
schemes, the solution of Euler equations is made possi- 
ble because of dissipative mechanisms introduced in the 
numerical methods. 

To conclude, although built up on empirical basis, the 
general EoS here proposed, shows the correct behaviour, 
even in the presence of dissipative shocks. 
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